Model Dependence of Lateral Distribution Functions of 
High Energy Cosmic Ray Air Showers 
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The influence of high and low energy hadronic models on lateral distribution functions of cosmic ray air 
showers for Auger energies is explored. A large variety of presently used high and low energy hadron interaction 
models are analysed and the resulting lateral distribution functions are compared. We show that the slope 
depends on both the high and low energy hadronic model used. The models are confronted with available 
hadron-nucleus data from accelerator experiments. 



I. INTRODUCTION 

Microscopic calculations of cosmic ray air showers at ener- 
gies > 10 18 eV are of crucial importance for the reconstruc- 
tion of the properties of the primary cosmic ray. Ground-array 
type experiments like AGASA and Auger (working in non- 
hybrid mode) determine the energy of the primary cosmic ray 
by measuring the particle density at some distance from the 
shower axis. Muons as ground particles are interesting be- 
cause they carry information about the composition of the pri- 
mary cosmic ray; showers induced by heavy nuclei tend to 
give more muons than showers induced by protons or pho- 
tons. 

The interpretation of the observed signals depends strongly 
on the quality of the air shower model and simulation. Careful 
modelling of the shower development and properties is there- 
fore of utmost importance. Many papers have restricted their 
focus to the dependence of air shower properties on the high 
energy hadronic model fliSiHli. This has been thought 
to be the main uncertainty for the determination of the shower 
maximum in the longitudinal profile. Models that predict a 
higher cross section or a higher multiplicity yield shower max- 
ima at smaller slant depths, i.e., higher in the atmosphere (see 
e.g. 0). 

In this paper we analyse the model dependence for lateral 
distribution functions (LDF) for muons, electrons/positrons 
and photons. For these observables, the low energy hadronic 
model becomes important, especially for the LDF of muons, 
because these are decay products of low energy charged 
mesons 0,0]. Also the tails of LDFs of electrons/positrons 
and photons are influenced by the low interaction energy re- 
gion as discussed in |9]. 



GHEISHA 0, G-FLUKAd (the notation G- denotes the 
older version of FLUKA as found in Geant3.21), and UrQMD 
1.2.1 1171 Il8l u% . Electromagnetic showering is done by the 
EGS4-codel2ai 

UrQMD is a microscopic transport approach based on the 
covariant propagation of constituent quarks and di-quarks ac- 
companied by mesonic and baryonic degrees of freedom. 
The leading hadrons of the fragmenting strings contain the 
valence-quarks of the original excited hadron. The elemen- 
tary hadronic interactions are modelled according to measured 
cross sections and angular distributions. If the cross sections 
are not experimentally known, detailed balance is employed 
in the energy range of resonances. The partial and total decay 
widths are taken from the Particle Data Group 1211 . 

Since UrQMD does not have an intrinsic cross section cal- 
culation, it is implemented in the following way. For colli- 
sions with e.g. nitrogen one assumes a disk with radius 5 fm, 
corresponding to a geometric cross section of 750 mb. The 
projectile is then propagated over a distance which is cho- 
sen randomly according to the mean free path, and is placed 
randomly on this disk. The fact that the cross section was 
assumed too large is automatically compensated by a certain 
number of non-interactions. The fraction of interactions times 
750 mb gives the inelastic cross section. In a certain way, one 
can picture this as nature does it: The particle passes through 
the air, and whenever the transverse distance to a nucleus is 
too large, it continues to propagate without interaction. 



III. RESULTS 



A. Modification of lateral distribution functions 



II. AIR SHOWER MODELLING 

For the present investigation of air shower models the 
SENECA framework is used. For an introduction and details 
about the SENECA model, the reader is referred to lIlolfTlll . 
In this approach, the intermediate energy region is calcu- 
lated by solving a set of transport equations numerically. The 
high energy models that we compare are QGSJET01 IU2I1 and 
SIBYLL2.1 dill; the choices for the low energy model are 



Lateral distribution functions are composed of three par- 
ticle types, namely electrons/positrons, photons and muons. 
The former two are initiated by the decays of high and low 
energetic 7r°s, whereas the muons are decay products of low 
energetic charged mesons. As pointed out in |9], high en- 
ergy pion decays determine the electromagnetic LDF near the 
core. However, at large distances from the core, the contribu- 
tion from low energy pions gains importance. Muons are in 
general produced by charged mesons that decay at intermedi- 
ate altitudes (~ 5 km). 
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Figure 1: Differences due to high energy interaction models for 
1 x 10 eV proton induced showers. Top panel: Ratio of LDFs of 
the model combinations (SIBYLL+GEISHA)/(QGSJET+GEISHA). 
Bottom panel: the longitudinal profiles (electrons/positrons and 
muons) of the same combinations. 



Let us start with a comparison of the LDFs for different 
high energy hadronic interaction models. We compare ver- 
tical proton induced 1 x 10 19 eV showers (averaged), ob- 
served at 860 g/cm 2 altitude. The kinetic energy cutoffs are 
1 MeV for electrons/photons, and 50 MeV for muons. As pre- 
viously shown with CORSIKA simulations in [5], the present 
study also gives fewer muons and more electrons for SIB YLL 
than for QGSJET at small distances. At large distances both 
the electromagnetic and the muonic component are smaller in 
SIB YLL, as shown in Fig.^Ctop). This is confirmed in the 
average longitudinal profile from both models, as shown in 
Fig- HI (bottom): QGSJET develops higher in the atmosphere. 
The significant difference in muons can be traced back to 
the higher average multiplicity of QGSJET. This difference is 
about 25% at 100 GeV and reaches a factor of two at 10 20 eV 
|0. At low energies this is due to different implementations 
of parton distribution functions (PDFs), QGSJET implements 
a flatter choice, even though contradictory to A flater PDF has 
a higher parton density a intermediate xf (Feynman-x: lon- 
gitudinal momentum fraction) or energy, and hi gher density 
at lower if ot higher energy. HERA data (seeTHH for 
discussion). This leads to a somewhat higher multiplicity at 
low energies (ca. 100 GeV) but helps to control particle mul- 
tiplicities up to collider energies. For higher energies parti- 
cle production increases again strongly since QGSJET misses 
higher order corrections. Current efforts focus on resolv- 
ing this problem by introducing enhanced Pomeron diagrams 




10000 



Figure 2: Dependence of LDFs for QGSJET as high energy model, 
and GHEISHA, G-FLUKA, and UrQMD as low energy models. 



|24|. SIB YLL mimics this effect by the implementation of an 
ad-hoc energy dependent p± cutoff for hard processes (3- 
Since the main contribution of particle production comes from 
hard processes, an increase of p± with energy leads to a mod- 
erate rise of the multiplicity. 

In the following, vertical 5 x 10 19 eV proton induced 
showers are analysed. Fig. [2] shows the ratio of LDFs for 
the following combination of hadronic models. QGSJET 
and GHEISHA (QG), QGSJET and G-FLUKA (QF), and 
QGSJET and UrQMD (QU). The functions are plotted as a ra- 
tio to QGSJET+UrQMD. Whereas the QGSJET+G-FLUKA 
and QGSJET+UrQMD models show rather similar LDFs, re- 
sulting in a ratio close to unity, QGSJET+GHEISHA has a 
notably flatter LDF The effect is most significant for muons, 
somewhat less for electrons and photons. This can be under- 
stood by the fact that most of the last collisions for muons 
producing hadrons are in the low-energy regime. The same 
is true for electrons and photons at larger distances from the 
shower-axis. The muon ratio reaches almost a factor of 2 at 
3000 meters. 

It is especially interesting to study the experimental re- 
sponse of the different components. At the Pierre Auger Ob- 
servatory, the density is measured in vertical equivalent muon 
units (VEM). A vertical minimum ionising muon counts as 
one particle. Since the detectors are water tanks, muons pen- 
etrate much deeper than electrons and photons, and are ex- 
pected to give a more pronounced signal. We calculate the 
detector response in terms of Cherenkov photons. A detector 
simulation based on the WTANK program 1 25 ] has been done 
for the setup of an Auger water tank, by counting the number 
of Cherenkov photons arriving at the positions of the photo- 
multipliers. Tables with the results have been made for differ- 
ent incident angles, energies and particle types and are then 
applied to the simulations. The results for an average vertical 
5 x 10 19 eV shower are shown in Fig. [3] The model combi- 
nation used here is QGSJETOl+UrQMD. One sees that the 
curves for muons in number density and VEM units almost 
agree, whereas photons and electrons give less contribution 
per particle. Photons are still dominant due to their high num- 
ber. These results are in good agreement with calculations 
presented in the Auger Design Report 1 26]. For completeness 
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Figure 3: Lateral distribution functions for different particle species 
in number densities and vertical equivalent muon units for an incom- 
ing proton of 5 x 10 19 eV. 
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Figure 4: LDF ratios of 5 x 10 eV proton induced vertical show- 
ers for QGS JET0 1 +GHEISH A and SIBYLL+G-FLUKA, the model 
combinations with the largest differences. The VEM of the total sig- 
nal is also shown. 



we show that protons and neutrons have a negligible contribu- 
tion to the total signal. 

The ratio of QGSJET+GHEISHA to SIB YLL+ G-FLUKA, 
together with the VEM yields, are depicted in Fig. |4] These 
model combinations give the most different results. The effect 
on the slope of the VEM-LDF is even bigger because at small 
distances photons dominate the signal, whereas at large dis- 
tances the muons give the most contribution. The sensitivity 
of Auger to models is somewhat higher than for experiments 
using thin plastic scintillators. The main reason is that wa- 
ter tanks are more sensitive to muons, since their mean path 
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Figure 5: The LDFs in VEM units, scaled to the fitted function for 
QGSJETOl+UrQMD. 



Model 


C 


V 


rn [m] 


a 


QGS JET0 1 +GHEISHA 


3.17xl0 10 


3.05 


986 


0.87 


QGSJETOl+UrQMD 


4.00 xlO 10 


3.10 


1540 


1.63 


QGS JET0 1 +G-FLUKA 


3.13xl0 10 


3.05 


1126 


1.34 


SIBYLL+GHEISHA 


4.88 xlO 10 


3.15 


1036 


0.96 


SIB YLL+UrQMD 


4.88 xlO 10 


3.15 


1136 


1.30 


SIBYLL+G-FLUKA 


4. 15 xlO 10 


3.11 


985 


1.26 



Table I: The fit parameters for function Q for all model combina- 
tions. 



length through the material is much longer. Further, muons 
become more dominant at large core distances, due to the flat- 
ter LDF. These two facts are important when choosing the 
density at a given distance as energy estimator. However, as 
shown below, GHEISHA as low energy hadronic model can 
be excluded from accelerator data and the differences between 
the remaining two low energy models are not so large. With 
a well measured LDF and a cross-calibration with the fluores- 
cence method, Auger might be able to discriminate between 
high energy hadronic models. 

The LDFs for the different models have been fitted to the 
following function 



S(r) = Cr-^l + (r/r ) 2 y 



(1) 



with C, 77, ro, a being fit parameters and r being the dis- 
tance to the shower axis in meter. This function is inspired 
from the one used by AGASA [27] . The resulting parameters 
are shown in Table |I] Fig. |5] shows the LDFs in VEM units 
for all model combinations, scaled by function Q, using the 
parameters for QGSJET+UrQMD. One sees how at small dis- 
tances, the high energy model determines the shape; at larger 
distances deviations result from differences in the low energy 
hadronic models. 



B. Comparisons of models with data 

In this section the behaviour of the different hadronic mod- 
els and their influence on the LDFs is analysed. Here we focus 
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Figure 6: Multiplicities for pions of a 7r +Ne reaction at 10.5 GeV/c. 
Data are taken from 11291. 



on data available for hadron-nucleus collisions in the relevant 
energy domain. Light nuclei are chosen which are somewhat 
comparable to air. Fig. [6] addresses the mean 4-7T multiplic- 
ity of pions for 7r~+Ne reactions at 10.5 GeV/c 1 28] . One 
clearly observes that GHEISHA produces too many pions, 
while UrQMD and G-FLUKA are in line with the measured 
data. 

To explore the longitudinal momentum distribution, we 
confront the models with rapidity spectra for p+Be reactions 
at 14.6 GeV/c taken at the AGS. The fact that GHEISHA 
shows a strange peak in the rapidity spectra is well-known 
13011 . In addition, GHEISHA yields too high multiplicities for 
pions, and G-FLUKA overestimates the longitudinal momen- 
tum loss for protons. The stopping of the initial proton is well 
described by UrQMD and GHEISHA. 

At higher energies, we can add QGSJET and SIBYLL to 
the comparison. These models are commonly employed at 
energies larger than £b cam ^100 GeV. For example, COR- 
SIKA uses 80 GeV as transition value. Fig. [8] shows that 
QGSJET overestimates and SIBYLL underestimates the neg- 
atively charged hadron multiplicities in p+Mg reactions at 100 
and 200 GeV/c. The three other model results are in line with 
the data (NA35 jH). 
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Figure 7: Rapidity spectra of final state pions and 
reactions at 14.6 GeV/c. Data are taken from E802 
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Figure 8: Multiplicity of negatively charged particles for a p+Mg 
reaction at 200 GeV/c. Data are taken from NA35 l3lll . 
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Figure 9: Rapidity distributions for the reaction p+Ar — > charged hadrons. Data are taken from l32ll . 




Figure 10: Rapidity distributions for the reaction p+Ar — > negatively charged hadrons. Data are taken from l32ll . 
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Figure 11: Rapidity distributions for the reaction p+S — > negatively charged hadrons. Data are taken from NA35 l33ll . 
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Figure 12: Rapidity distributions for the reaction p+S — » net protons. Data are taken from NA35 l33ll . 
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Figure 13: Mean multiplicities of negatively charged hadrons for the 
discussed set of low energy models. 
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Figure 14: Mean multiplicities of negatively charged hadrons for the 
discussed set of high energy models. UrQMD results are given as a 
comparison. 



Figs.|9]andrro|show rapidity spectra of (negatively) charged 
hadrons for p+Ar reactions 1 32] . Again, QGSJET and 
SIBYLL do deviate from the data. GHEISHA still produces 
the peak, which is located such that the total multiplicity is 
well described. G-FLUKA and UrQMD give good descrip- 
tions of the data. A similar behaviour is obtained for p+S 
collisions at 200 GeV/c (NA35 data |'330), as demonstrated in 
Figs. and [21 Note that the analysis has been performed 
with the experimental trigger demanding a minimum of 5 



charged particles being produced in the reaction. 

As above, the longitudinal momentum loss of the protons 
(stopping) is best described by G-FLUKA and UrQMD. 

An overall comparison is given in Figs. 1131 and !14l where 
the mean multiplicity of negatively charged hadrons is plotted 
as a function of the incident kinetic energy. Below 100 GeV 
GHEISHA produces more particles per collision; this leads to 
the higher muon densities at large distances. In the same way 
the multiplicity of produced 7r°s changes the shape of the LDF 
of electrons and photons, only at larger distances, where this 
energy region becomes more dominant. UrQMD is slightly 
lower in multiplicity than G-FLUKA. The fact that the lat- 
ter results in smaller signals at large distances might be due 
to differences in the total cross section, which is neglected 
in this analysis. The overestimation of particle production by 
QGSJET in the transition region above 100 GeV is significant. 
Since QGSJET is the most used model in air shower simula- 
tions, these deviations do deserve further studies, because the 
energy region above 100 GeV is certainly important for muon 
production. 



IV. CONCLUSIONS 

• We have demonstrated that the LDFs for high energy 
air showers depend on the high and also on the low en- 
ergy hadronic models. The tails are strongly influenced 
by the low energy model, which therefore modifies the 
slope of the LDFs. 

• GHEISHA, which is commonly used as low energy 
hadronic model, reproduces less accurately available 
data of hadron induced interactions. G-FLUKA and 
UrQMD describe the data better. 

• In the transition region around 100 GeV, QGSJET pre- 
dicts too high particle multiplicities, whereas SIBYLL 
is somewhat low. The present findings suggest to in- 
crease the transition energy. UrQMD is tested up to 
RHIC energies (21 TeV), and could be employed up to 
this value. 
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